package fem2.material;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import math2.MathUtilities;
import math2.MutableDouble;
import fem2.Debugger;
import fem2.Element;
import fem2.MeshPart;
import fem2.Variable;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.State;
import fem2.enu.VariableName;

/**
 * Drucker-Prager plasticity model without kinematic hardening for:<br>
 * + Plane strain and 3d<br>
 * + Non-associative flow rule<br>
 * + Associative isotropic hardening<br>
 * + Implicit backward Euler & return mapping
 * 
 * This model follow the formulation in Souze-Neto book but with different
 * incremental setting.
 * 
 * @author hbui
 * 
 */
public class DruckerPragerMM extends InelasticMaterialModel {

	private double eta, etabar, xi;
	private FluencyCriteria hardeningLaw;

	/**
	 * 
	 * @param E
	 *            elastic modulus
	 * @param nu
	 *            Poisson ratio
	 * @param thickness
	 *            thickness
	 * @param rho
	 *            density
	 * @param ss
	 *            Stress state
	 * @param eta
	 *            friction coeeficient
	 * @param etabar
	 *            non-associative flow friction coeeficient
	 * @param xi
	 *            dilatancy factor
	 * @param hardeningLaw
	 *            isotropic hardening law
	 */
	public DruckerPragerMM(double E, double nu, double thickness, double rho, State ss, double eta,
			double etabar, double xi, FluencyCriteria hardeningLaw) {
		super(E, nu, thickness, rho, ss);

		if (getState() == State.PLANE_STRESS) {
			throw new Error(
					"Because of the complexity, this material model does not support plane stress version. Please use DruckerPragerPsMM instead");
		}

		this.eta = eta;
		this.etabar = etabar;
		this.xi = xi;
		this.hardeningLaw = hardeningLaw;
	}

	@Override
	public void Initialize(Element e) {
		super.Initialize(e);

		int dim = e.getMeshPart().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		int ntrs1 = (dim == 2 ? ntrs + 1 : ntrs);

		ArrayList<Variable> vars = getVariables(e);
		for (Variable v : vars) {
			v.registerScalarVariable(VariableName.ALPHA, true); // strain-like
																// variable
			v.registerVectorVariable(VariableName.PLASTIC_STRAIN, ntrs1, true);
			/*
			 * save the variable to create history
			 */
			v.save();
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		StructuralElement se = (StructuralElement) e;
		MeshPart mp = e.getMeshPart();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		int ntrs1 = (dim == 2 ? ntrs + 1 : ntrs);
		int ng = e.getNumIntegrationPoints();

		IVector eetrial = new ArrayVector(ntrs1);
		IVector eedtrial = new ArrayVector(ntrs1);
		IVector strial = new ArrayVector(ntrs1);

		if (dim == 2) {

			/*************************************************************
			 * Drucker-Prager plane strain
			 *************************************************************/

			for (int k = 0; k < ng; k++) {
				double[] xi = e.getIntegrationPoint(k);
				Variable var = getVariables(e).get(k);
				double alphaOld = var.getScalarVariable(VariableName.ALPHA, -1);
				MutableDouble alpha = var.getScalarVariable(VariableName.ALPHA);
				IVector epOld = var.getVectorVariable(VariableName.PLASTIC_STRAIN, -1);
				IVector ep = var.getVectorVariable(VariableName.PLASTIC_STRAIN);
				IVector epsilon = var.getVectorVariable(VariableName.STRAIN);
				IMatrix C = getConstitutiveMatrix(e, k);
				IVector sigma = se.getStress(k);

				BLAM.copy(se.computeStrain(xi), epsilon);

				// Debugger.watch("epsilon = ", epsilon);

				/*
				 * compute trial elastic strain
				 */
				eetrial.set(0, epsilon.get(0) - epOld.get(0));
				eetrial.set(1, epsilon.get(1) - epOld.get(1));
				eetrial.set(2, epsilon.get(2) - epOld.get(2));
				eetrial.set(3, -epOld.get(3)); // for plane strain
				// epsilon(3)==0

				/*
				 * compute trial deviatoric elastic strain
				 */
				double eevtrial = (eetrial.get(0) + eetrial.get(1) + eetrial.get(3));
				eedtrial.set(0, eetrial.get(0) - eevtrial / 3);
				eedtrial.set(1, eetrial.get(1) - eevtrial / 3);
				eedtrial.set(2, eetrial.get(2));
				eedtrial.set(3, eetrial.get(3) - eevtrial / 3);

				/*
				 * compute trial deviatoric stress.
				 */
				strial.set(0, 2 * mu * eedtrial.get(0));
				strial.set(1, 2 * mu * eedtrial.get(1));
				strial.set(2, mu * eedtrial.get(2));
				strial.set(3, 2 * mu * eedtrial.get(3));

				/*
				 * compute trial pressure
				 */
				double ptrial = kappa * eevtrial;
				// Debugger.watch("ptrial = ", ptrial);

				/*
				 * check yield condition
				 */
				double J2 = (strial.get(0) * strial.get(0) + strial.get(1) * strial.get(1) + 2
						* strial.get(2) * strial.get(2) + strial.get(3) * strial.get(3)) / 2;
				double sqrtJ2 = Math.sqrt(J2);
				double c = hardeningLaw.computeFluencyValue(alphaOld);

				double ftrial = sqrtJ2 + this.eta * ptrial - this.xi * c;

				if (ftrial / c < MathUtilities.tol12) {

					/*
					 * update internal variables as trial state
					 */
					alpha.setValue(alphaOld);
					BLAM.copy(epOld, ep);

					/*
					 * update stress
					 */
					sigma.set(0, strial.get(0) + ptrial);
					sigma.set(1, strial.get(1) + ptrial);
					sigma.set(2, strial.get(2));

					/*
					 * update tangent operator
					 */
					BLAM.copy(this.De, C);
				} else {

					/*
					 * return to smooth portion of the cone: compute deltaLambda
					 */
					double dlambda = computeDeltaLambdaSmooth(sqrtJ2, ptrial, alphaOld);

					double smoothValidation = sqrtJ2 - mu * dlambda;

					// if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
					// Debugger.watch("sqrtJ2 -mu*dlambda= ", smoothValidation);
					// }

					// Debugger.watch("dlambda = ", dlambda);

					/*
					 * check validity condition
					 */

					if (smoothValidation >= 0) {

						/*
						 * returning point lie in smooth portion of the cone
						 */

						if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
							Debugger.watch("return to smooth cone portion");
						}

						/*
						 * update alpha
						 */
						alpha.setValue(alphaOld + this.xi * dlambda);

						/*
						 * compute N/sqrt(2) to update plastic strain directly.
						 * Enable this when updating plastic strain directly
						 * below
						 */
						// IVector N2 = (IVector) strial.clone();
						// BLAM.scale(0.5 / Math.sqrt(J2), N2);

						/*
						 * update deviatoric stress
						 */
						BLAM.scale(1 - mu * dlambda / sqrtJ2, strial);

						/*
						 * update pressure
						 */
						ptrial -= kappa * etabar * dlambda;

						/*
						 * update stress
						 */
						sigma.set(0, strial.get(0) + ptrial);
						sigma.set(1, strial.get(1) + ptrial);
						sigma.set(2, strial.get(2));

						/*
						 * update elastic strain to update plastic strain
						 * indirectly
						 */
						eetrial.set(0, strial.get(0) / (2 * mu) + ptrial / (3 * kappa));
						eetrial.set(1, strial.get(1) / (2 * mu) + ptrial / (3 * kappa));
						eetrial.set(2, 2 * strial.get(2) / (2 * mu)); //
						eetrial.set(3, strial.get(3) / (2 * mu) + ptrial / (3 * kappa));

						/*
						 * update plastic strain
						 */
						ep.set(0, epsilon.get(0) - eetrial.get(0));
						ep.set(1, epsilon.get(1) - eetrial.get(1));
						ep.set(2, epsilon.get(2) - eetrial.get(2));
						ep.set(3, -eetrial.get(3));

						/*
						 * update plastic strain directly
						 */
						// BLAM.copy(epOld, ep);
						// ep.add(0, dlambda * (N2.get(0) + etabar / 3));
						// ep.add(1, dlambda * (N2.get(1) + etabar / 3));
						// ep.add(2, dlambda * 2 * N2.get(2));
						// ep.add(3, dlambda * (N2.get(3) + etabar / 3));

						/*
						 * update tangent operator
						 */
						// IVector D = (IVector) eedtrial.clone();
						IVector D = eedtrial;
						D.set(2, eedtrial.get(2) / 2);// convert back to tensor
														// form
						double norm_eedtrial = Math.sqrt(D.get(0) * D.get(0) + D.get(1) * D.get(1)
								+ 2 * D.get(2) * D.get(2) + D.get(3) * D.get(3));

						if (norm_eedtrial > MathUtilities.tol6) {
							BLAM.scale(1.0 / norm_eedtrial, D);
						} else {
							BLAM.scale(0.0, D);
						}

						double H = hardeningLaw.computeFluencyDerivative(alpha.getValue());
						double A = 1.0 / (mu + kappa * eta * etabar + this.xi * this.xi * H);
						double aFactor = 2 * mu * (1 - dlambda / (sqrt_2 * norm_eedtrial));
						double bFactor = 2 * mu * (dlambda / (sqrt_2 * norm_eedtrial) - mu * A);
						double cFactor = -sqrt_2 * mu * A * kappa;
						double dFactor = kappa * (1 - kappa * eta * etabar * A);

						// IMatrix DbyD = new Array2DMatrix(ntrs, ntrs);
						// LinAlgUtilities.multiply(1.0, D, D, DbyD);
						//
						// IMatrix DbyI = new Array2DMatrix(ntrs, ntrs);
						// DbyI.set(0, 0, D.get(0));
						// DbyI.set(0, 1, D.get(0));
						// DbyI.set(1, 0, D.get(1));
						// DbyI.set(1, 1, D.get(1));
						// DbyI.set(2, 0, D.get(2));
						// DbyI.set(2, 1, D.get(2));
						//
						// IMatrix IbyD = new Array2DMatrix(ntrs, ntrs);
						// IbyD.set(0, 0, D.get(0));
						// IbyD.set(1, 0, D.get(0));
						// IbyD.set(0, 1, D.get(1));
						// IbyD.set(1, 1, D.get(1));
						// IbyD.set(0, 2, D.get(2));
						// IbyD.set(1, 2, D.get(2));
						//
						//
						// BLAM.zero(C);
						// BLAM.add(aFactor, LinAlgUtilities.Rank4Identity2d,
						// C);
						// BLAM.add(-aFactor / 3, LinAlgUtilities.OneByOne2d,
						// C);
						// BLAM.add(bFactor, DbyD, C);
						// BLAM.add(cFactor * eta, DbyI, C);
						// BLAM.add(cFactor * etabar, IbyD, C);
						// BLAM.add(dFactor, LinAlgUtilities.OneByOne2d, C);

						for (int i = 0; i < ntrs; i++) {
							for (int j = 0; j < ntrs; j++) {
								// C.set(i,
								// j,
								// aFactor
								// * LinAlgUtilities.FOID[i][j]
								// + bFactor
								// * D.get(i)
								// * D.get(j)
								// + cFactor
								// * (eta * D.get(i) * LinAlgUtilities.SOID[j] +
								// etabar
								// * D.get(j) * LinAlgUtilities.SOID[i])
								// + (dFactor - aFactor / 3) *
								// LinAlgUtilities.SOID[i]
								// * LinAlgUtilities.SOID[j]);
								C.set(i,
										j,
										aFactor
												* LinAlgUtilities.Idev2d.get(i, j)
												+ bFactor
												* D.get(i)
												* D.get(j)
												+ cFactor
												* (eta * D.get(i) * LinAlgUtilities.SOID[j] + etabar
														* D.get(j) * LinAlgUtilities.SOID[i])
												+ dFactor * LinAlgUtilities.SOID[i]
												* LinAlgUtilities.SOID[j]);
							}
						}

					} else {

						/*
						 * return to apex
						 */
						if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
							Debugger.watch("return to cone apex");
						}

						/*
						 * return to the apex: compute depv
						 */
						double depv = computeDepvApex(ptrial, alphaOld);
						dlambda = depv / etabar;

						/*
						 * update alpha
						 */
						alpha.setValue(alphaOld + this.xi * dlambda);

						/*
						 * update pressure
						 */
						ptrial -= kappa * etabar * dlambda;

						/*
						 * update sigma
						 */
						sigma.set(0, ptrial);
						sigma.set(1, ptrial);
						sigma.set(2, 0);

						// /*
						// * update elastic strain
						// */
						// eetrial.set(0, ptrial / (3 * kappa));
						// eetrial.set(1, ptrial / (3 * kappa));
						// eetrial.set(2, 0);
						// eetrial.set(3, ptrial / (3 * kappa));
						//
						// /*
						// * update plastic strain indirectly
						// */
						// ep.set(0, epsilon.get(0) - eetrial.get(0));
						// ep.set(1, epsilon.get(1) - eetrial.get(1));
						// ep.set(2, epsilon.get(2) - eetrial.get(2));
						// ep.set(3, -eetrial.get(3));

						/*
						 * update plastic strain directly
						 */
						BLAM.copy(epOld, ep);
						ep.add(0, depv * etabar / 3);
						ep.add(1, depv * etabar / 3);
						ep.add(3, depv * etabar / 3);

						/*
						 * update tangent operator
						 */
						double H = hardeningLaw.computeFluencyDerivative(alpha.getValue());
						double factor = kappa
								* (1 - kappa / (kappa + this.xi * this.xi * H / (eta * etabar)));

						BLAM.copy(LinAlgUtilities.OneByOne2d, C);
						BLAM.scale(factor, C);

						// for (int i = 0; i < ntrs; i++) {
						// for (int j = 0; j < ntrs; j++) {
						// C.set(i, j, factor * LinAlgUtilities.SOID[i]
						// * LinAlgUtilities.SOID[j]);
						// }
						// }

					}
				}

				if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
					// Debugger.watch("sigma = ", sigma);
					if (e.getId() == 1) {
						Debugger.watch("uHatLocal = ", e.getMeshPart().getUHatLocal());
					}
				}

			}
		}

		/*************************************************************
		 * Drucker-Prager 3d
		 *************************************************************/

		else if (dim == 3) {

			for (int k = 0; k < ng; k++) {
				double[] xi = e.getIntegrationPoint(k);
				Variable var = getVariables(e).get(k);
				double alphaOld = var.getScalarVariable(VariableName.ALPHA, -1);
				MutableDouble alpha = var.getScalarVariable(VariableName.ALPHA);
				IVector epOld = var.getVectorVariable(VariableName.PLASTIC_STRAIN, -1);
				IVector ep = var.getVectorVariable(VariableName.PLASTIC_STRAIN);
				IVector epsilon = var.getVectorVariable(VariableName.STRAIN);
				IMatrix C = getConstitutiveMatrix(e, k);
				IVector sigma = se.getStress(k);

				BLAM.copy(se.computeStrain(xi), epsilon);

				// Debugger.watch("epsilon = ", epsilon);

				/*
				 * compute trial elastic strain
				 */
				BLAM.copy(epsilon, eetrial);
				BLAM.add(-1.0, epOld, eetrial);

				/*
				 * compute trial deviatoric elastic strain
				 */
				double eevtrial = (eetrial.get(0) + eetrial.get(1) + eetrial.get(2));
				eedtrial.set(0, eetrial.get(0) - eevtrial / 3);
				eedtrial.set(1, eetrial.get(1) - eevtrial / 3);
				eedtrial.set(2, eetrial.get(2) - eevtrial / 3);
				eedtrial.set(3, eetrial.get(3));
				eedtrial.set(4, eetrial.get(4));
				eedtrial.set(5, eetrial.get(5));

				/*
				 * compute trial deviatoric stress.
				 */
				strial.set(0, 2 * mu * eedtrial.get(0));
				strial.set(1, 2 * mu * eedtrial.get(1));
				strial.set(2, 2 * mu * eedtrial.get(2));
				strial.set(3, mu * eedtrial.get(3));
				strial.set(4, mu * eedtrial.get(4));
				strial.set(5, mu * eedtrial.get(5));

				/*
				 * compute trial pressure
				 */
				double ptrial = kappa * eevtrial;
				// Debugger.watch("ptrial = ", ptrial);

				/*
				 * check yield condition
				 */
				double J2 = (strial.get(0) * strial.get(0) + strial.get(1) * strial.get(1)
						+ strial.get(2) * strial.get(2) + 2 * (strial.get(3) * strial.get(3)
						+ strial.get(4) * strial.get(4) + strial.get(5) * strial.get(5))) / 2;
				double sqrtJ2 = Math.sqrt(J2);
				double c = hardeningLaw.computeFluencyValue(alphaOld);

				double ftrial = sqrtJ2 + this.eta * ptrial - this.xi * c;

				if (ftrial / c < MathUtilities.tol12) {

					/*
					 * update internal variables as trial state
					 */
					alpha.setValue(alphaOld);
					BLAM.copy(epOld, ep);

					/*
					 * update stress
					 */
					sigma.set(0, strial.get(0) + ptrial);
					sigma.set(1, strial.get(1) + ptrial);
					sigma.set(2, strial.get(2) + ptrial);
					sigma.set(3, strial.get(3));
					sigma.set(4, strial.get(4));
					sigma.set(5, strial.get(5));

					/*
					 * update tangent operator
					 */
					BLAM.copy(this.De, C);
				} else {

					/*
					 * return to smooth portion of the cone: compute deltaLambda
					 */
					double dlambda = computeDeltaLambdaSmooth(sqrtJ2, ptrial, alphaOld);

					double smoothValidation = sqrtJ2 - mu * dlambda;

					// if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
					// Debugger.watch("sqrtJ2 -mu*dlambda= ", smoothValidation);
					// }

					// Debugger.watch("dlambda = ", dlambda);

					/*
					 * check validity condition
					 */

					if (smoothValidation >= 0) {

						/*
						 * returning point lie in smooth portion of the cone
						 */

						if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
							Debugger.watch("return to smooth cone portion");
						}

						/*
						 * update alpha
						 */
						alpha.setValue(alphaOld + this.xi * dlambda);

						/*
						 * compute N/sqrt(2) to update plastic strain directly.
						 * Enable this when updating plastic strain directly
						 * below
						 */
						// IVector N2 = (IVector) strial.clone();
						// BLAM.scale(0.5 / Math.sqrt(J2), N2);

						/*
						 * update deviatoric stress
						 */
						BLAM.scale(1 - mu * dlambda / sqrtJ2, strial);

						/*
						 * update pressure
						 */
						ptrial -= kappa * etabar * dlambda;

						/*
						 * update stress
						 */
						sigma.set(0, strial.get(0) + ptrial);
						sigma.set(1, strial.get(1) + ptrial);
						sigma.set(2, strial.get(2) + ptrial);
						sigma.set(3, strial.get(3));
						sigma.set(4, strial.get(4));
						sigma.set(5, strial.get(5));

						/*
						 * update elastic strain to update plastic strain
						 * indirectly
						 */
						eetrial.set(0, strial.get(0) / (2 * mu) + ptrial / (3 * kappa));
						eetrial.set(1, strial.get(1) / (2 * mu) + ptrial / (3 * kappa));
						eetrial.set(2, strial.get(2) / (2 * mu) + ptrial / (3 * kappa));
						eetrial.set(3, 2 * strial.get(3) / (2 * mu));
						eetrial.set(4, 2 * strial.get(4) / (2 * mu));
						eetrial.set(5, 2 * strial.get(5) / (2 * mu));

						/*
						 * update plastic strain
						 */
						BLAM.copy(epsilon, ep);
						BLAM.add(-1.0, eetrial, ep);

						/*
						 * update plastic strain directly
						 */
						// BLAM.copy(epOld, ep);
						// ep.add(0, dlambda * (N2.get(0) + etabar / 3));
						// ep.add(1, dlambda * (N2.get(1) + etabar / 3));
						// ep.add(2, dlambda * (N2.get(2) + etabar / 3));
						// ep.add(3, dlambda * 2 * N2.get(3));
						// ep.add(4, dlambda * 2 * N2.get(4));
						// ep.add(5, dlambda * 2 * N2.get(5));

						/*
						 * update tangent operator
						 */
						// IVector D = (IVector) eedtrial.clone();
						IVector D = eedtrial;
						D.set(3, eedtrial.get(3) / 2);// convert back to
														// tensor form
						D.set(4, eedtrial.get(4) / 2);
						D.set(5, eedtrial.get(5) / 2);
						double norm_eedtrial = Math.sqrt(Math.pow(D.get(0), 2)
								+ Math.pow(D.get(1), 2)
								+ Math.pow(D.get(2), 2)
								+ 2
								* (Math.pow(D.get(3), 2) + Math.pow(D.get(4), 2) + Math.pow(
										D.get(5), 2)));

						if (norm_eedtrial > MathUtilities.tol6) {
							BLAM.scale(1.0 / norm_eedtrial, D);
						} else {
							BLAM.scale(0.0, D);
						}

						double H = hardeningLaw.computeFluencyDerivative(alpha.getValue());
						double A = 1.0 / (mu + kappa * eta * etabar + this.xi * this.xi * H);
						double aFactor = 2 * mu * (1 - dlambda / (sqrt_2 * norm_eedtrial));
						double bFactor = 2 * mu * (dlambda / (sqrt_2 * norm_eedtrial) - mu * A);
						double cFactor = -sqrt_2 * mu * A * kappa;
						double dFactor = kappa * (1 - kappa * eta * etabar * A);

						for (int i = 0; i < ntrs; i++) {
							for (int j = 0; j < ntrs; j++) {
								// C.set(i,
								// j,
								// aFactor
								// * LinAlgUtilities.Rank4Identity3d.get(i, j)
								// + bFactor
								// * D.get(i)
								// * D.get(j)
								// + cFactor
								// * (eta * D.get(i) * LinAlgUtilities.SOID3d[j]
								// + etabar
								// * D.get(j) * LinAlgUtilities.SOID3d[i])
								// + (dFactor - aFactor / 3)
								// * LinAlgUtilities.SOID3d[i]
								// * LinAlgUtilities.SOID3d[j]);
								C.set(i,
										j,
										aFactor
												* LinAlgUtilities.Idev3d.get(i, j)
												+ bFactor
												* D.get(i)
												* D.get(j)
												+ cFactor
												* (eta * D.get(i) * LinAlgUtilities.SOID3d[j] + etabar
														* D.get(j) * LinAlgUtilities.SOID3d[i])
												+ dFactor * LinAlgUtilities.SOID3d[i]
												* LinAlgUtilities.SOID3d[j]);
							}
						}

					} else {

						/*
						 * return to apex
						 */
						if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
							Debugger.watch("return to cone apex");
						}

						/*
						 * return to the apex: compute depv
						 */
						double depv = computeDepvApex(ptrial, alphaOld);
						dlambda = depv / etabar;

						/*
						 * update alpha
						 */
						alpha.setValue(alphaOld + this.xi * dlambda);

						/*
						 * update pressure
						 */
						ptrial -= kappa * etabar * dlambda;

						/*
						 * update sigma
						 */
						sigma.set(0, ptrial);
						sigma.set(1, ptrial);
						sigma.set(2, ptrial);
						sigma.set(3, 0.0);
						sigma.set(4, 0.0);
						sigma.set(5, 0.0);

						// /*
						// * update elastic strain
						// */
						// eetrial.set(0, ptrial / (3 * kappa));
						// eetrial.set(1, ptrial / (3 * kappa));
						// eetrial.set(2, ptrial / (3 * kappa));
						// eetrial.set(3, 0);
						// eetrial.set(4, 0);
						// eetrial.set(5, 0);
						//
						// /*
						// * update plastic strain indirectly
						// */
						// BLAM.copy(epsilon,ep);
						// BLAM.add(-1.0, eetrial, ep);

						/*
						 * update plastic strain directly
						 */
						BLAM.copy(epOld, ep);
						ep.add(0, depv * etabar / 3);
						ep.add(1, depv * etabar / 3);
						ep.add(2, depv * etabar / 3);

						/*
						 * update tangent operator
						 */
						double H = hardeningLaw.computeFluencyDerivative(alpha.getValue());
						double factor = kappa
								* (1 - kappa / (kappa + this.xi * this.xi * H / (eta * etabar)));

						BLAM.copy(LinAlgUtilities.OneByOne3d, C);
						BLAM.scale(factor, C);

						// for (int i = 0; i < ntrs; i++) {
						// for (int j = 0; j < ntrs; j++) {
						// C.set(i, j, factor * LinAlgUtilities.SOID3d[i]
						// * LinAlgUtilities.SOID3d[j]);
						// }
						// }

					}
				}

			}
		}

	}

	/**
	 * compute delta lambda for return map to smooth portion of cone
	 * 
	 * @param sqrtJ2
	 * @param alpha
	 * @return
	 */
	private double computeDeltaLambdaSmooth(double sqrtJ2, double ptrial, double alpha) {
		double dlambda = 0.0;
		double a = alpha;
		double p = ptrial;
		int cnt = 0;
		while (++cnt < MAX_NT) {
			double c = hardeningLaw.computeFluencyValue(a);
			double g = sqrtJ2 - mu * dlambda + eta * p - xi * c;

			if (Math.abs(g / c) < MathUtilities.tol8) {
				break;
			}

			double H = hardeningLaw.computeFluencyDerivative(a);
			double dg = mu + kappa * eta * etabar + xi * xi * H;
			dlambda += g / dg;
			a = alpha + xi * dlambda;
			p = ptrial - kappa * etabar * dlambda;

			// Debugger.watch("sqrtJ2 =", sqrtJ2);
			// Debugger.watch("ptrial =", ptrial);
			// Debugger.watch("alpha =", alpha);
			// Debugger.watch("a = ", a);
		}
		if (cnt == MAX_NT) {
			throw new Error("dlambda computation does not converge in " + MAX_NT + " loops");
		}
		return dlambda;
	}

	/**
	 * compute delta epv for return map to apex
	 * 
	 * @param normxi
	 * @param alpha
	 * @return
	 */
	private double computeDepvApex(double ptrial, double alpha) {
		double depv = 0.0;
		double a = alpha;
		double p = ptrial;
		int cnt = 0;
		while (cnt++ < MAX_NT) {
			double c = hardeningLaw.computeFluencyValue(a);
			double g = (xi / eta) * c - p;

			if (Math.abs(g / c) < MathUtilities.tol8) {
				break;
			}

			double H = hardeningLaw.computeFluencyDerivative(a);
			double dg = xi * xi * H / (eta * etabar) + kappa;
			depv -= g / dg;
			a = alpha + (xi / etabar) * depv;
			p = ptrial - kappa * depv;
		}
		if (cnt == MAX_NT) {
			throw new Error("dlambda computation does not converge in " + MAX_NT + " loops");
		}
		return depv;
	}
}
